rm(list =ls())
options(scipen=999)
gc()
packages <-c("tidyverse","estimatr","plm","stargazer",
             "fastDummies","ICCbin","ihs","readstata13","xtable")

new.packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new.packages)) install.packages(new.packages)

lapply(packages, require, character.only = TRUE)
rm(packages, new.packages)

setwd("PUT YOUR DIRECTORY HERE")

## insheet data
data <- read.csv("./Datasets/officers_longitudinal_reshaped_ub.csv")

## scale outcome
data$ub_staff <- scale(data$ub_staff)
data$commander_ <- scale(data$commander_)

## run models
fe_model_year <- plm(ub_staff ~ commander_, data = data, index = c("year"), model = "within")
fe_model_year_locality <- plm(ub_staff ~ commander_, data = data, index = c("year", "locality"), model = "within", effect = "twoways")


stargazer(fe_model_year, fe_model_year_locality, 
          dep.var.labels = c("UB Personell", "UB Personell"),
          style = "qje",
          covariate.labels = c("Secret Police \n Officers"),
          star.char = c("*", "**", "***"),
          star.cutoffs = c(0.1, 0.05, 0.01),
          omit = c("Constant"),
          omit.stat = c("rsq", "f", "ser", "adj.rsq"), 
          omit.table.layout = "n",
          add.lines = list(c("Time FEs", "Yes", "YES"),
                           c("Unit FEs", "NO", "Yes")),
          out = "PUT YOUR FILEPATH HERE")
